Accurate time— domain gravitational waveforms for extreme-mass-ratio binaries 
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The accuracy of time-domain solutions of the inhomogeneous Teukolsky equation is improved 
significantly. Comparing energy fluxes in gravitational waves with highly accurate frequency- 
domain results for circular equatorial orbits in Schwarzschild and Kerr, we find agreement to within 
1% or better, which may be even further improved. We apply our method to orbits for which 
frequency-domain calculations have a relative disadvantage, specifically high-eccentricity (elliptical 
and parabolic) "zoom-whirl" orbits, and find the energy fluxes, waveforms, and characteristic strain 
in gravitational waves. 
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I. INTRODUCTION 

A stellar mass compact object (a black hole or a neu- 
tron star) moves in the spacetimc of a supermassive black 
hole in an (accelerated) orbit, emitting low- frequency 
gravitational waves in the good sensitivity band of space 
borne detectors such as LISA. In this Paper we model a 
likely scenario in galactic nuclei, specifically the capture 
of compact stellar-size objects by a supermassive black 
hole, like those that exist in the center of many galax- 
ies. Because of the high mass ratio, the compact object 
may be approximated by a point particle, and the orbital 
evolution modeled by perturbation theory. 

The emission of gravitational waves by a point-particle 
falling into a black hole has been modeled by several re- 
searchers using a frequency-domain (FD) decomposition 
of the master equation (see, e.g., [1, 2] and references 
cited therein) [3] . The FD approach allows for very accu- 
rate determination of the energy flux and the waveforms 
for wide classes of orbits that include orbits character- 
ized by low values of the eccentricity. Generic orbits and 
especially high eccentricity orbits, although in principle 
amenable to a Fourier decomposition and a FD construc- 
tion of the waveforms [2, 4], require the summation of 
many terms in the Fourier series, which limits the accu- 
racy and increases the computation time [5]. This weak- 
ness of the FD approach is not shared with the time- 
domain (TD) approach [6, 8, 9] that we present here. 
Indeed, wave generation is typically one of the strengths 
of TD approaches. We expect a TD approach to be more 
efficient and yield more accurate results (for waveforms) 
especially in the context of highly elliptic and parabolic 
orbits compared with FD calculations. TD calculations 
can also serve as an independent test of FD results. 

Calculations in the TD have only recently been in- 
troduced for the cases of equatorial circular particle or- 
bits [8] and also for equatorial elliptic orbits and inclined 
circular orbits [6] around a rapidly rotating black hole. 
These calculations, however, suffer from low accuracy, as 
was manifest in deviations by as much as 10-20% in the 



determination of the energy flux when compared with 
the highly accurate FD results. Indeed, it is this crudity 
of TD calculations that was the main reason why so few 
TD calculations have been done. Here, we report on a 
significant improvement in the accuracy of TD calcula- 
tions, that brings them to within 1% agreement with the 
FD counterparts or better. Then, we apply our approach 
to a special class of orbits, namely "zoom-whirl" orbits 
of both elliptic [10] and parabolic type. As was recently 
pointed out in [11], adiabatic waveforms that ignore the 
gravitational wave backreaction on the orbit may yield 
so-called "kludge" waveforms, that are surprisingly ac- 
curate for much of the parameter space of extreme-mass- 
ratio binary motion. We adopt here a similar approach, 
and neglect the effects of radiation reaction, hoping to 
address them elsewhere. The purpose of this Paper is to 
serve as a proof-of- concept that TD calculations of grav- 
itational waveforms can be accurate — specifically where 
FD calculations have weaknesses — and hence be of prac- 
tical importance for waveform generation. This Paper is 
not intended to serve as a reference for tables of results, 
or to exhaustively cover the relevant parameter space. 
The latter — in addition to additional improvements to 
the code and model — will be presented elsewhere. 

This paper is organized as follows. In Section II we 
briefly review the TD approach. In Section III we demon- 
strate the greatly improved accuracy of our approach 
by reproducing results obtained with FD calculations to 
within 1%. Next, in Section IV, we present new results 
for accurate TD waveforms for equatorial elliptic and 
parabolic, zoom-whirl orbits. Finally, in Section V we 
discuss our results. 



II. TIME-DOMAIN CALCULATIONS 

The Teukolsky equation, that governs the perturba- 
tions of Kerr geometry in Boyer-Lindquist coordinates 
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with a matter-source term, is given by [12] 
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All the symbols used above are defined in [12]. By choos- 
ing their values in Boyer-Lindquist coordinates and ex- 
pressing Eq. (1) explicitly we obtain an explicit but quite 
complicated form for the source term. We treat the delta 
functions that appear in (2) by using a "regularization" in 
which we substitute the delta function with a very nar- 
row Gaussian function (represented only by a few grid 
points), 



5(x - x(t)) 



1 



2tt (7 



exp 



2o 2 



(3) 



for small a. Wc handle the deltas in the radial and 
angular direction through the Gaussian approximation, 
whereas the 5((j) — <j>(t)) function is handled analytically. 
The methodology for numerically evolving the homoge- 
neous part of (1) has been presented in detail [13]. We 
take a similar approach, except that we include the source 
term which is essential for modeling extreme-mass-ratio 
binaries. More details can be found in [8]. 

Note that we take the initial data for the fields to be 
zero. Since we are considering the Tcukolsky equation 
with a source, this choice corresponds to the particle "ap- 
pearing suddenly," which generates an artificial burst of 
spurious radiation, that decays fast. We handle this diffi- 
culty by evolving to late times and computing fluxes only 
after the system has "settled down." We used the highly 
accurate public domain integrator "Geod" [14] — which 
we call from within our Tcukolsky code — to update the 
particle's position at every time step. 

Finally, to obtain the power ("total energy flux") due 
to gravitational wave emission from the binary, we find 
the Teukolsky function tp at a set radial location in the 
grid far from the horizon, and compute [15] 
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The characteristic strain in gravitational waves is then 
estimated by 
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where z is the cosmological redshift, Dl is the luminosity 
distance, and | dE e / df e \ is the total energy per unit fre- 
quency carried away by the gravitational waves (in the 
"emitter's" reference frame, hence the subscripts "e"). 



III. TESTING THE CODE'S ACCURACY 

The numerical implementation of this work is in the 
form of a set of 2+1 dimensional linear PDE's that uses 
the Lax-Wendroff evolution scheme. This approach, in- 
cluding convergence and stability, is described in [8] . 

We test the accuracy of our code by performing stan- 
dard convergence tests, and we find the code to converge 
with second order. We also tested the code by comparing 
our results for the energy flux in gravitational waves with 
those obtained with FD calculations. Specifically, we 
compare circular and equatorial orbits for Schwarzschild 
and Kerr black holes, where the FD fluxes are known to 
high accuracy. 

As discussed above, earlier TD calculations suffered 
from crude accuracy, with energy fluxes typically devi- 
ating from their FD counterparts by 10-20%. Analysis 
of the reasons for that crudity allowed us to attribute it 
primarily to the relatively coarse grids used in [6, 8], and 
to the relatively small size of the computational domain, 
which required the evaluation of the field only a small 
number of wavelengths from the source, in addition to 
unoptimized modeling of the point particle. Moreover, 
there exists a typographical sign error in the expression 
for Cm,n in Ref. [16] which was used in [6, 8] for com- 
puting the source-term. This wrong sign introduces a 
significant error in the strong-field region of the black 
hole space-time, but has little impact when orbits with 
large radii are considered. Finally, it should also be noted 
that the extra normalization factor as used before (see 
Eq. 20 in [8]) is actually inappropriate in the context of 
our computation. In our current code we have fixed the 
sign error as mentioned above and removed this improper 
normalization. Also, owing to certain optimization and 
performance improvements to the code, we can now have 
much larger grid sizes — and therefore extract the wave- 
forms and compute the fluxes at a distance much farther 
than in [6, 8] — and also have much higher resolution. 

In Fig. 1A we show the relative error in the deter- 
mination of the energy flux in gravitational waves as a 
function of the distance at which wave extraction is done, 
for a particle of mass \xjM = 0.01 in a circular equatorial 
orbit of radius 5ri sco about a Kerr black hole with spin 
parameter a/M — 0.9. In practice, we extract the wave- 
forms and compute the fluxes at or beyond 500M. Here, 
and similarly throughout this paper, the "detector" is 
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FIG. 1: The relative error in the energy flux in gravitational 
waves for a particle of mass (J./M = 0.01 in a circular equa- 
torial orbit of radius 5ri sco about a Kerr black hole with spin 
parameter a/M — 0.9. Upper panel (A): As a function of 
the distance at which wave extraction is done. The errors 
are calculated with a value corresponding to wave extraction 
at infinity, that we obtain using Richardson's extrapolations. 
Here, N = 5. Lower panel (B): As a function of the number 
of points used to sample the Gaussian N. The errors are cal- 
culated with the FD value. Wave extraction is done at 500M. 
Data for both panels are calculated at our chosen grid density. 



positioned on the equatorial plane of the black hole, al- 
though there are no restrictions in positioning it at any 
inclination angle with respect to the equatorial plane. 
Our typical grid resolutions are 0.025M (radial) and 0.05 
(angular) . Modeling the particle by a Gaussian distribu- 
tion requires optimization of the number of grid points 
N that are used to sample the Gaussian: For a given grid 
resolution, too many grid points cause the physical width 
of the Gaussian to spread too much resulting in large de- 
viations from the point-particle approximation. On the 
other hand, too few grid points cause the Gaussian to be 
"under-sampled," introducing unwanted errors into the 
computation. Through extensive numerical experimen- 
tation, we have found the optimal number of grid points 
across the gaussian for our chosen grid resolution to be 
about five. Figure IB shows the relative errors in the 
determination of the flux extracted at 500M as a func- 
tion of N, for the same physical system. For the modes 
to = 2, 3 that we report here, we find that N — 5 is 
optimal. This value for N may change for other values 
for to. We also experimented in modeling the particle 
with non-Gaussian profiles, e.g., high-order polynomials 
with compact support, and a "discrete-delta" approach 
werein one starts out by defining a step function on a 
numerical grid and takes finite-differenced derivatives to 
obtain an approximation to the delta function and its 
derivatives. Done optimally, the results one obtains from 
these various implementations are consistent. Further 
details about the "discrete-delta" approach will be dis- 
cussed elsewhere. We consider this as one of the strengths 



TABLE I: Comparison of gravitational wave total energy 
fluxes detected at infinity using the FD solution, with the 
results measured numerically at r — 500M while evolving 
the Teukolsky equation in the TD. The particle of mass 
fi/M = 0.01 is in a circular equatorial orbit of radius r — n sco 
about a Schwarzschild black hole. 
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7.31 x 10~ 8 


7.37 x 10~ 8 
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1.45 x 10~ 8 


1.46 x 10~ 8 
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3.57 x 10~ 9 


3.60 x 10~ 9 



TABLE II: Same as Table I, for a particle of mass n/M = 0.01 
in a circular equatorial orbit of radius 5.0n SC o about a Kerr 
black hole with spin parameter a/M = 0.9. 
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2.19 x 10~ 9 


2.21 x 10~ 9 
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2.20 x 10~ 10 
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2.68 x 10~ n 


2.71 x 10~ n 



of our approach: our results are quite robust. 

Implementing these refinements in the code, we consid- 
erably improve the accuracy of our TD results compared 
with the previous standard. Sample comparisons with 
FD are described in tables I and II. Our results agree 
with FD results to about 1% or better. We anticipate 
that this agreement can be further improved, by an ad- 
ditional increase in resolution (possibly using a multi-grid 
approach of extrapolation to infinite resolution) , extract- 
ing the field farther out (again, with a possible extrap- 
olation to infinity), and by more experimentation with 
various models for the point particle. 
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FIG. 2: Equatorial zoom-whirl orbits used in the text. Left 
panel (A): An elliptic orbit about a Kerr black hole with spin 
a/M = 0.5, with p = 5.0 M and e = 0.5. This orbit in- 
cludes a "whirl" part that has about four full quasi-circular 
cycles, before it enters the "zoom" phase. Right panel (B): A 
parabolic orbit about a Kerr black hole with spin a/M = 0.5, 
with p = 5.828427 M and e = 1. This orbit includes a "whirl" 
part that has about eight full quasi-circular cycles, before it 
enters the "zoom out" phase. 
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FIG. 3: Waveform from the elliptic orbit. We display Re(?/)) 
and lm(tp) as functions of time. Note the "initial burst" which 
should be ignored. The smooth, slow oscillations from the 
"whirl" part of the orbit are clearly visible. The waveform 
drops off very rapidly in the "zoom" part of the orbit. Panel 
A: Real part of the dominant mode (m = 2). Panel B: Imagi- 
nary part of the same. Panel C: Real part of the mode m — 3. 
Panel D: Imaginary part of the same. 



IV. EQUATORIAL ZOOM-WHIRL ORBITS 

"Zoom-whirl" orbits have been studied in detail [10]. 
The basic idea is that for certain elliptic orbits (close 
to the separatrix between stable and unstable orbits) the 
amount of time spent close to the periastron is very large, 
so that the particle traces a quasi-circular path close to 
the periastron before traveling back to the apastron. For 
high-eccentricity orbits, the particle appears to "zoom 
in" towards its periastron location, spend a number of 
circular revolutions ("whirls"), and then "zoom out" to- 
wards the apastron. These orbits are very interesting, 
especially for the case of rapidly spinning black holes, be- 
cause of the large number of "whirls" they perform close 
to the hole. Similarly, there are parabolic orbits that 
exhibit similar "whirl" and "zoom out" characteristics. 

First, we plot sample equatorial zoom- whirl orbits 
about a Kerr black hole with spin a/M — 0.5 in Fig. 2A 
(elliptic orbit with semi-latus rectum p = 5.0 M and ec- 
centricity e = 0.5) and in Fig. 2B (parabolic orbit with 
p = 5.828427 M and e = 1). 

We next show the waveforms and energy fluxes from a 
particle moving in the orbits depicted in Figs. 2A (in 
Figs. 3-4, henceforth the "elliptic orbit") and 2B (in 
Figs. 5-7, hereafter the "parabolic orbit"). The wave- 
forms presented are the real and imaginary parts of the 
m = 2, 3 modes of the Teukolsky function. All waveforms 
are extracted at 500M (i.e., these are the waveforms for 
an observer at that distance.) As mentioned above, the 
initial burst of radiation is spurious, as it corresponds 
to the particle "appearing out of nowhere." The inter- 
esting features of these waveforms are as follows: The 



LU 




t/M 



FIG. 4: Total energy flux from the elliptic orbit for the dom- 
inant mode m = 2 (solid curve) and the m = 3 mode (dashed 
curve). Again, the "initial burst" should be ignored. The al- 
most steady flux from the "whirl" part of the orbit is clearly 
visible. 
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FIG. 5: Waveform from the parabolic orbit for the m — 2 
mode (upper panel). We display Re(i/>) (solid curve) and 
Im(ip) (dashed) as functions of time. The spurious "initial 
burst" is plotted in the lower panel. The smooth, slow oscil- 
lations from the "whirl" part of the orbit are clearly visible. 
The waveform drops off very rapidly in the "zoom out" part 
of the orbit. 



emitted waves become very prominent as the particle 
approaches the periastron, and fall off rapidly when the 
particle moves away from the hole. Note that the part of 
the waveform from the "whirl" seems to be similar to one 
coming from an exact circular orbit of the same angular 
velocity. In fact, we checked this quantitatively, and in- 
deed, the waveform of the "whirl" part exactly matches 
with the waveform that arises from a corresponding circu- 
lar orbit. This happens both for the elliptic and parabolic 
case. 

It is clear from the figures that the physically signif- 
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FIG. 6: Same as Fig. 5, for the mode m = 3. 
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FIG. 7: Total energy flux from the parabolic orbit for the 
dominant mode m = 2 (solid curve) and the m = 3 mode 
(dashed curve). The initial burst was excised. The almost 
steady flux from the "whirl" part of the orbit is clearly visible. 

icant parts of zoom whirl orbits, as expected, are the 
"whirls" because the emitted radiation from that part is 
much stronger. This is interesting, because this implies 
that with some care, one may be able to estimate radia- 
tive fluxes etc. for such orbits using known circular orbit 
results. Notice that although the parabolic orbit is sym- 
metric (about the x-axis, i.e., the "zoom in" phase of he 
orbit is symmetrical to the "zoom out" phase) , the wave- 
forms and the energy flux are not symmetrical. Specifi- 
cally, notice the oscillations in (the lower panel of) Fig. 7, 
that are visible only during the "zoom out" phase. The 
reason for this asymmetry is that during the "zoom out" 
phase of the orbit, the waves include also the waves emit- 
ted during the "zoom in" and the "whirl" phases, that 
backscatter off the curvature of spacetime. Indeed, the 
(real part of the) angular frequency of the oscillations in 
Fig. 7 is found numerically to equal ~ 0.4 M^ 1 , which 
is very close to twice the orbital angular frequency of 




FIG. 8: Characteristic strain (bold curve) in gravitational 
waves h c for a IMq black hole moving on the parabolic orbit 
around a 1O 6 M0 black hole with a/M = 0.5, at a distance of 
lGpc, shown on a standard LISA noise curve with signal-to- 
noise ratio of 1 [17]. 

0.183 M -1 . (Note that the decay of the "initial burst" is 
characterized by oscillations with the quasinormal mode 
frequency.) Finally, we show in Fig. 8 the characteristic 
strain in gravitational waves h c for the parabolic orbit, 
for the dominant mode m — 2. 

Notice that while for circular orbits (both 
Schwarzschild and Kerr) the energy flux in gravita- 
tional waves in the dominant mode m = 2 is greater 
than that in the next mode (m = 3) by an order of 
magnitude, for the zoom-whirl orbits we have studied it 
is greater by only a factor of ~ 2. This difference may 
have implications on the number of modes that needs to 
be computed. 

The new results we are presenting here for zoom-whirl 
orbits can be expected to be accurate within a 1% or 
better. However, as mentioned before, we expect to be 
able to reduce this error through higher resolution, better 
particle modeling, and related ideas. 



V. CONCLUSIONS 

A variety of particle orbits in the context of extreme- 
mass-ratio binaries can be effectively studied using Kerr 
black hole perturbation theory in the TD. We demon- 
strate the potential accuracy of the TD approach by get- 
ting very good agreement (within 1% or better) in com- 
parisons with FD results for radiated energy fluxes for 
a particle in a circular orbit. This agreement may be 
further improved. 

We also obtained gravitational "waveforms" arising 
from particles in a special class of orbits, namely zoom- 
whirl orbits of both elliptic and parabolic type. This is 
the first study of Kerr parabolic orbits. We observed that 
the "whirl" part of these orbits dominates the emitted ra- 
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diation and is very similar to the pattern of radiation aris- 
ing from a circular orbit of the appropriate parameters. 
We also present the characteristic strain in gravitational 
waves for a parabolic zoom- whirl orbit. We ignored the 
effects of radiation reaction in this work, with the hope 
of addressing those in detail elsewhere. 

The TD Teukolsky approach to extreme-mass-ratio in- 
spiral is maturing, and no longer suffers from the poor 
accuracy problems that used to plague it. One may ex- 
pect interesting and useful results from it where the FD 
approach has disadvantages, in addition to being a useful 



independent check on FD results. 
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